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The transition from Townsend to glow discharge is investigated numerically in one space dimen- 
sion in full parameter space within the classical model: with electrons and positive ions drifting 
in the local electric field, impact ionization by electrons (a process), secondary electron emission 
from the cathode (7 process) and space charge effects. We also perform a systematic analytical 
small current expansion about the Townsend limit up to third order in the current that fits our 
numerical data very well. Depending on the two determining parameters 7 and system size pd, the 
transition from Townsend to glow discharge can show the textbook subcritical behavior, but for 
smaller values of pd, we also find supercritical or some unexpected intermediate "mixed" behavior. 
Our work shows the same qualitative dependence of i7 = U{I,pd) for fixed 7 as the old experiments 
by Pokrovskaya-Soboleva and Klyarfeld. Furthermore, the analysis lays the basis for understanding 
the complex spatio-temporal patterns in short planar barrier discharge systems. 
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I. INTRODUCTION 

Space charge effects in many cases are the first non- 
linear effects in gas discharges with increasing current. 
They are known to induce the avalanche to streamer 
transition in transient discharges as well as the transi- 
tion from Townsend to normal and further to abnormal 
glow in stationary discharges. Generically, nonlinear cou- 
plings in non-equilibrium systems lead to the formation 
of spontaneous spatio-temporal patterns. The current 
constriction in the normal glow discharge as well as the 
longitudinal striations of a long positive column of a glow 
discharge [1-3] fall into this class of phenomena. 

Recently, the amazing variety of spatio-temporal pat- 
terns formed mainly in the transversal direction of a short 
dc driven system consisting of a gas discharge layer and 
a semiconductor layer sandwiched between two copla- 
nar electrodes has drawn considerable attention [4-12]. 
These patterns are due to the nonlinear gas discharge 
being coupled to the linearly responding semiconduc- 
tor. In particular, a negative differential conductivity of 
the gas discharge in some region of the current-voltage- 
characteristics is expected [13-19] to play a significant 
role in the spontaneous formation of patterns, quite like 
in nonlinear semiconductor devices [20]. Due to its ge- 
ometry, modeling the system [10-12] as one-dimensional 
is a very good approximation, as long as this symmetry 
is not spontaneously broken by the intrinsic dynamics. 
So as a first step of any investigation, the behavior and 
the resulting current- voltage-characteristics of the purely 
one-dimensional gas discharge system have to be under- 
stood. 

An investigation of the system [11] along the lines of 
the textbook [21] shows that the pattern formation oc- 
curs at the space charge driven transition from Townsend 



to glow discharge. The gas dicharge layer is rather short, 
more precisely, the product pd of gas pressure p times 
electrode distance d is small. This raises the question 
of the Townsend to glow transition for small pd. How- 
ever, despite a history of more than 70 years, we are not 
aware of any thorough and complete study of this classi- 
cal problem. Therefore, our aim in the present paper is 
to develop a consistent picture of the Townsend to glow 
transition in one dimension from analytical and numeri- 
cal investigations, in particular, for short systems. 

Many authors focus on quite long discharges that have 
a clearly pronounced subcritical characteristics, i.e., for 
fixed large pd and growing total current /, the voltage 
first decreases from the Townsend limit towards the nor- 
mal glow regime, then it increases again in the abnormal 
glow regime until heating effects become important and 
the voltage again decreases towards the arc discharge. 
We will not consider this last thermally driven transi- 
tion at high currents. The initial decrease of voltage 
from Townsend discharge towards normal glow creates 
a regime of negative differential conductivity, and some 
authors [19] believe that negative differential conductiv- 
ity is generic for this system. 

However, already in the early 1940'ies, e.g., in the ex- 
tensive review by Druyvesteyn and Penning [22], it was 
suggested that this subcritical behavior might not be the 
only possible one, but that also a monotonic increase 
of voltage with current was possible. Such a behavior 
we will call supercritical, in line with modern bifurca- 
tion theory. There are early experimental papers by 
Pokrovskaya-Soboleva and Klyarfeld [23] and McClure 
[24] that clearly indicate a supercritical transition for 
small values of pd in hydrogen and deuterium in com- 
bination with metal electrodes. Later data by the same 
authors [25] is reproduced in Raizer's textbook [21], how- 
ever, only for rather long systems with subcritical char- 
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acteristics. 

Theoretical insight into the question of bifurcation be- 
havior can be gained by analytical or numerical inves- 
tigation of the appropriate model. The classical model 
contains the drift of charged particles in the local field, 
the a-process of impact ionization in the bulk of the gas, 
the 7-process of secondary electron emission from the 
cathode, and space charge effects. 

Numerical calculations date back to the 50'ics [26], 
the first numerical evaluations using an "electronic com- 
puter" can be found in the early 60'ies in [27,28]. In par- 
ticular, in the work of Ward [28], current- voltage charac- 
teristics with or without a region of negative differential 
conductivity can be found for different values ofpd. How- 
ever, computing power at the time was quite restricted 
and hence only a few current- voltage-characteristics were 
calculated. The work does not seem to have been ex- 
tended significantly lateron. We will take up the issue in 
Section IV. 

Analytical efforts were constrained to small current ex- 
pansions about the Townsend limit. The old German 
textbook of Engel and Steenbeck [13] contains an ele- 
gant argument that the initial increase or decrease of the 
characteristics from the Townsend limit depends on the 
sign of a"{ET) where a{E) is the effective impact ioniza- 
tion coefficient as a function of the electric field E, and 
" denotes the second derivative evaluated at Townsend's 
breakdown field Et- We recall this argument in Sec- 
tion III.B. The book [13] also gives an explicit expres- 
sion for the coefficient oc a" {Ex) in the expansion 
U{I) = Ut + C2/^, however, without derivation or refer- 
ence. Exactly the same statements can be found more 
than 60 years later in Raizer's much read textbook [21]. 
Kolobov and Fiala [29] assume that a" = marks the 
point where negative differential conductivity disappears. 
A similar small current expansion of the voltage about 
the Townsend limit has recently been performed in [19], 
but with a different result — here the leading correction 
is found to be linear in the current rather than quadratic. 
None of the two results has been compared to numerical 
solutions. In the present paper, we will present yet an- 
other result for the small current expansion and evaluate 
it to higher orders. Our derivation is a systematic ex- 
pansion and in very good agreement with our numerical 
results. 

In general, our aim in the present paper is a con- 
sistent theoretical investigation of the simple classical 
model of these discharges treated by so many authors 
[13,18,19,21,22,26-30]. The exploration of the full pa- 
rameter space is possible, because the current-voltage 
characteristics in appropriate dimensionless units de- 
pends essentially only on two parameters: the secondary 
emission coefficent 7 and the dimensionless system size 
L (X pd. 

Of course, various extensions of the model can be con- 
sidered: particle diffusion, attachement, nonlinear parti- 
cle mobilities, a field-dependent secondary emission rate 
or nonlocal ionization rates. However, e.g., Boeuf [31] 



has argued that for the transition from normal to abnor- 
mal glow, nonlocal terms in the impact ionization reac- 
tion should be taken into account through hybrid numer- 
ical models [32], while in the subnormal regime between 
Townsend and normal glow, a local fluid model is con- 
sidered sufficient [29] . This supports the strategy to first 
seek a full understanding of the predictions of the clas- 
sical model as a corner stone and starting point for any 
further work. 

In the present work, we perform a systematic analyt- 
ical expansion of the voltage about the Townsend limit 
up to 0{P), recovering the qualitative features of the 
solution from [13,21]: in particular, we find that a linear 
term in current / indeed is missing, and that the coef- 
ficient C2 indeed is proportional to a"{ET), but with a 
different proportionality constant. In fact, our coefficient 
C2 depends strongly on the secondary emission coefficient 
7 — it varies by almost three orders of magnitude for 7 
between 10~^ to 10~^ — , while the expression given in 
[13,21] does not depend on 7 at all. We also evaluate the 
next order 0{I^). Our analytical result fits our numer- 
ical solutions very well within its range of validity. The 
stationary states of the pattern forming system [11] are 
within the range of validity of this expansion. 

Furthermore, we explore the current-voltage character- 
istics numerically beyond the range of the small current 
expansion in the full parameter space. We show that 
within the classical model, there is not only the familiar 
subcritical bifurcation from Townsend to glow discharge 
for large values of pd, but for sufficiently small values 
of pd, the bifurcation is supercritical, in agreement with 
the scenario suggested by Druyvesteyn and Penning [22]. 
Furthermore, for intermediate values of pd, there always 
exist completely unexpected "mixed" bifurcations. This 
surprising finding implies that the negative differential 
conductivity does not vanish when a"{ET) = in the 
Townsend limit, as most authors assume [21,29], but only 
for smaller values of pd. These statements are true for 
all relevant values of secondary emission 7. Our three- 
dimensional plots of the voltage as a function of dimen- 
sionless system size L oc pd and current / for a given 
gas-electrode combination are done in the same manner 
as the old experimental plots by Pokrovskaya-Soboleva 
and Klyarfeld [25]. 

The paper is organized as follows: in Section II, we 
recall the classical model and its parameters, perform di- 
mensional analysis, and reformulate the stationary one- 
dimensional problem as a boundary condition problem. 
In Section HI, we recall the Townsend limit and the clas- 
sical argument of Engel and Steenbeck on the qualitative 
dependence of the small current expansion on a". We 
then perform a new systematic small current expansion 
up to third order in the total current and determine 
the coefficients of the expansion explicitly. Section IV 
begins with our numerical strategy and a discussion of 
the parameters with their ranges. The parameter de- 
pendence of the current- voltage characteristics on system 
size L oc pd and secondary emission coefficient 7 is first 



2 



presented in the form of (/, U,pd)-j>\ots for fixed 7 as in 
[25]. We then present spatial plots of electron current 
and field, and compare our numerical results to our an- 
alytical small current expansion. Finally, we classify the 
bifurcation structure in the complete relevant parameter 
space. Section V contains a summary and an outlook 
onto the implications of this work for spatio-temporal 
pattern formation in barrier discharges. Two appendices 
contain the proof of the uniqueness of the solution of the 
boundary value problem and details of the small current 
expansion in order I^. 



is characterized by the single parameter s with typical 
values s = 1/2 or 1 depending on the type of gas. Our 
numerical results are for the most common value s = 1. 

The parameter 7 is the probability that a positive ion 
hitting the cathode leads to emission of a free electron 
into the gas. For a discharge of length d with the anode 
at X = and the cathode at X = d, the 7 process enters 
as boundary condition a.t X = d 

\Jeid,t)\ = J \J+{d,t)\ , (7) 

while ions are absent at the anode 



II. THE CLASSICAL MODEL 

A. Definition 

We investigate the classical model for glow discharges 
in simple non-attaching gases in a planar, quasi-one- 
dimensional geometry. The same model was previously 
investigated in. e.g., [13,21,22,26-30] as discussed in the 
introduction. The model consists of continuity equations 
for two charged species, namely electrons and positive 
ions with particle densities Ug and n+ 



dt Ue + dxJe = source , 
dt n+ + dxJ+ = source . 



(1) 
(2) 



Their space charges can modify the externally applied 
field E through the Poisson equation 



dxE = — (n+ - Tie) . 

So 



(3) 



In the simplest approximation, diffusion is neglected and 
particle current densities Jg and J+ are approximated by 
drift only 



-Ue HeE , J+ =n+ IJ,+ E 



(4) 



where the drift velocity here is assumed to be linearly 
dependent on the local field with mobilities <C /ie- 

Two ionization processes are taken into account: the 
a process of ionization by electron impact in the bulk of 
the gas, and the 7 process of electron emission by ion 
impact onto the cathode. In a local field approximation, 
the a process is modeled as a local source term in the 
continuity equations 



J+(0,i)=0. (8) 
The electric potential U between the electrodes is 

U{t) ^ $(0, t) - t) > , E{X, t) = -dx^ ■ (9) 

With this convention, the average electric field E is pos- 
itive. Equations (l)-(9) define the classical model. 

B. Reformulation and dimensional analysis 

For the further calculation, it is useful to note, that 
the continuity equations (1), (2) together with the Pois- 
son equation (3) in one dimension result in the spatial 
conservation of the total electric current 

eodtE + e{J+-Je) = J{t) , dxJ = 0, (10) 

and that the ion current density J+ (4) with the help of 
(3) can be completely expressed by Jg and E 



J+ = !^(-Je+'-^t^eEdxE) . 
Me V e / 



(11) 



By dimensional analysis, the independent dimension- 
less parameters of the model are identified. It is con- 
venient to introduce the following dimensionless times, 
lengths and fields 



X 



t 



to ' 



= , .(.,.) (12) 



no 



En 



where 



source = \Je\ a{\E\) , a{\E\) = Ap a 



Bp 



(5) 



where p is the pressure of the gas. (The mobilities then 
scale with inverse pressure fig = fie/p and = j2+/p.) 
In the classical Townsend approximation [21], the func- 
tion 



a(£)=e-(VI^I)^ 



(6) 



Xo = -^ , Eo=Bp, (13) 

— = HeEo = IJ-eB p , no = -^rr- = p . 

to eXo e 



The equations now take the form 

dr£ = j{T) -{1 + IJ.)je - nSd^S , 



(14) 
(15) 
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where je = (t£ = — cJe/[cno Xo/to] is the dimensionless 
conductive current carried by the electrons, 



J 



oc 



eno Xo/to p' 



and u 



U 



oc 



(16) 



are the dimensionless total current and potential, and 





fj, = — oc p 

Me 



and L= = Apd (17) 
Xo 



are the ratio of ion over electron mobility and the dimen- 
sionless length of the gas discharge layer. 

We here have also recalled the scaling properties with 
pressure p, such that the pressure dependent similarity 
laws easily can be identified in the dimensionless results 
below. 



C. The stationciry problem 



jeix) =j e 



(24) 



and by evaluating this solution with the boundary con- 
dition je{L) = j e~^i at L. The identity (23) also can 

be found in [13,21]. 

It follows immediately that for a bounded function 
with Q.{£) < 1 for all £ as in (6), the system size L 
needs to be larger than 



L> 



(25) 



to sustain a stationary self-sustained discharge. This is 
true for arbitrary currents j and space charge effects. 

The identity (23) also plays a prominent role in the 
small current expansion about the Townsend limit, as we 
will see now. 



For a given dimensionless total current j, mobility ra- 
tio /i, secondary emission coefficient 7, functional form 

a(£) as in (6) and dimensionless system length L, the 
stationary solutions of (14), (15) are determined by 



dxje = -a{£) je , 
n£dx£ = J - (1 + ^J■)je 



(18) 
(19) 



together with the boundary conditions (7), (8) that con- 
veniently are expressed by j as 



je(0)=j and je{L)=je 



with 



In 



1 + 7 



(20) 



(21) 



We assume that a{£) > and da/d\£\ > within the 
relevant range of fields £. We prove in Appendix A that 
this determines a unique solution for the two functions 
je{x) and £{x). Finally, the integrated field yields the 
potential 

u = [ £{x) dx , (22) 
and hence the current- voltage-characteristics u{j). 

D. A global conservation law 

a{£{x)) for all solutions {jp{x), £{x)) is related to Lj 
through the global conservation law 



a{£{x)) dx = . 



(23) 



This can be seen by formally integrating Eq. (18) with 
the boundary condition je{0) = j with the result 



III. ANALYTICAL SMALL CURRENT 
EXPANSION 



A. The Townsend limit 



The well-known Townsend limit can be understood as a 
consequence of (23): for currents j so small that dx£ ~ 
in (19), the electric field is constant £{x) = £t- Eq. (23) 
then reduces to the familiar "ignition condition" [21] 



a{£T)L = 



(e - 1) = 1 . (26) 



The Paschen curve relates the potential ut = £tL in the 
Townsend limit to the system size L through a{uT/L) = 
L^/L. In particular, for the form of Eq. (6), the Paschen 
curve is 

T 

UT{L,'y) 

while the field is 

£t{L,^) 



\n^l'{L/L^) ' 



(27) 



(28) 



In dimensionless form, ut and £t depend only on the 

secondary emission coefficient 7, system size L and the 
parameter s in (6). The Townsend field £t increases 
monotonically with decreasing system size L and diverges 
for L [ L^. The Paschen curve ut{L,'j) (27) has a min- 
imum at L = e^/" and diverges both for L I and 
for L — > 00. 
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B. The argument of Engel and Steenbeck 

In the old German textbook of Engel and Steenbeck 
[13], the following argument for an expansion about the 
Townsend limit can be found: write the electric field as 
the Townsend field £t plus a perturbation A (a;), and note 
that the potential is the integrated field: 

£{x) =£t + A(x) , u = ut+ [ A(x) dx . (29) 

Jo 

The local impact ionization coefficient can then be ex- 
panded about a{£T) as 



a{£{x)) = a{£T) + o'C^t) A(x) + 



— A^{x) + . . 



(30) 



For fixed system size L and parameter L^, the global 
constraint (23) relates difi'erent solutions £{x) to a{£T)L 
through 

a{£T)L = I a{£{x)) dx 
Jo 

= a{£T)L + a' (£t) i A{x) dx 
Jo 



+ 



a"i£T) 



I A^{x) dx + ... , 
Jo 



(31) 



where the expansion of a was used in the second step. 
This identity allows to express A{x) dx by the higher 

order terms A"(x) dx, n = 2,3, — Insertion of this 
expansion into the definition of u yields 



Ut 



a"{£T) 
2 a'{£T) 



[ A'^{x)dx + ... 
Jo 



(32) 



Since A^ is positive and since a is assumed to be an in- 
creasing function of £, the sign of the correction is deter- 
mined by the sign of a". This statement from [13] is re- 
called in the recent literature [21,29]. It should be noted 
that the estimate (32) is valid as long as | a^"^ / A'^dx | -C 
\a"jA'^dx\ foraUn>3. 



anode: £{L) « 0. This prescription yields no dependence 
on 7 at all, quite in contrast to our results below. The 
functional forms for £{x) should be compared with our 
systematic analytical results (33), (49) below (note that 
we reversed the order of anode and cathode), and with 
our numerically derived field profiles in Figs. 4 and 5. 
They do not justify the ansatze given above. 

Rather a consistent ansatz is chosen in [19], and the 
structure of their expansion in terms of e~^"' and is 
quite similar to ours below. However, these authors fail 
to incorporate the global conservation law (23), and get 
a correction already in linear order, in contrast to the 
rigorous result (32) above. 



C. A systematic expansion in small j 

We now perform a systematic, expansion in powers of 
j about the Townsend limit. In principle, this expansion 
can be extended to arbitrary order. We have evaluated 
it up to 0{j^). We write the field correction as a power 
series in j, namely A (a;) = j £i{x) + £2{x) + . . ., and 
use the same ansatz for the current je{x) 



£{x) =£t+j £i (x) + f £2{x) + ... , 

je{x) = j il{x) +f L2{x) + ... , 

and we introduce the short hand notation 

a = a{£T) , a = a{£T) , a" = a"{£T) 
for the Taylor expansion of 
a{£{x)) =a + a' {j£i{x) + f£2{x) + . . .) 

+ ^(jfi(x)+j2£:2(x) + ...)' 



(33) 
(34) 



(35) 



(36) 



Insertion of the ansatze (33), (34) into Eqs. (18) and or- 
dering in powers of j yields 



0{j^) : d^ii{x) 
0{f) : d^L2{x) 



-ii{x) a , (37) 
- L2{x) a — i\{x) a' £\{x) , ... (38) 



The question is now how to caknilate A'^{x)dx. In For Eq. (19), the same procedure gives 



[13], a result is quoted referring to a long calcidation 
without reference. The same result is given more than 
60 years later in [21] in Section 8.3 with a sketch of an 
argument and again without reference. The argument as- 
sumes that \J+ \ 3> \Je \ throughout the discharge volume. 
This assumption is in disagreement with the boundary 
condition (8). A somewhat difi'erent argument based on a 
constant space charge through the whole system is given 
in [23]. In Ref. [21], the electric field profile is assumed 
to be £{x) oc i/l — x/xq, while [23] it is assumed to be 
£{x) cx (1 — x/xq) where the length scale .xq depends on 
the current j. In both cases, the breakdown of the ap- 
proximation is determined from the field vanishing at the 



o{f) 

0{f) 



d^£T = , (39) 
^i£Td^£i ^ I - {I + ^l)il{x) , (40) 
^£Tdx£2 + tJ'£idx£i = -(1 + M)t2(a;), • • • (41) 



The boundary condition (20) at the anode (x = 0) yields 
ti(0) = l , t2(0) = , t3(0)=0 , ... (42) 

The boundary condition (20) at the cathode {x = L) 

most conveniently is evaluated with the help of the global 
conservation law (23). Taking into account that L-y is in- 
dependent of j, the expanded form reads 
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where the first equation (43) reproduces the ignition con- 10 _g _g _^ ^ _2 .| , 

dition (26). Finally, the potential u from (22) is ""O ''^ ""O ""O lOyiO 



u = UT{L,-f)+j £i{x)dx + j^ £2{x) dx 
Jo Jo 



+ f I Ssix) dx + ... 



L 



(47) 



The lowest order ut{L,^) reproduces the Paschen curve 
(27). Eq. (44) reveals immediately that the order in u 
has to be absent. For the order p in (47), the function 
£i{x) has to be calculated. First, 



ii{x) = e 



(48) 



is the solution of (37) and (42). b\{x) has to be inserted 
into (40) which now can be solved analytically up to a 
constant of integration. This constant is determined by 
(44). The result is 

For the contribution in order p to the potential, the cal- 
culation of £i is sufficient since with the help of (45): 



/ £2[x) dx = — - — - / t-t (x) dx = — - — 7 — 
Jo 2a' Jo ^ 2a' a 



a" F(7,M) 



with the function 



^(7, M) = ^ + (1 + m) (2 - - 2e-^- - L^e-^-) 



+ (i + m) 



2,1-6-2^^ (1-6-^^)2 



. (50) 



The function is plotted in Fig. 1. Within the interesting 
parameter regime, it depends strongly on 7 and invisibly 
on fi. Here wc use the parameter range for 7 suggested by 
[21] and the maximal mobility ratio n = n+j ^Jte = 0.0095 
is reached for the lightest molecules, namely hydrogen. 



FIG. 1. Plot of -^(7, /i) as a function of 7 in a double- logarithmic 
plot. The dependence on ^, for realistic values Q < n < 0.0095 is 
too weak to be visible in the plot. However, F(7,/j) varies over 
almost 4 orders of magnitude as a function of 7. 

The small current expansion of the current-voltage- 
characteristics is in this approximation 



U = Ut 



I^J 2 a' {a£T) 



0{3^). (51) 



The range of validity of this expansion can be easily es- 
timated by inserting (49) into (33): the correction to the 
field due to the current should not exceed half of the 
Townsend field, so 



£t 



7<i £t /x^ ^^2) 



2max|£:i(x)| 2 £i{L) 



In view of the very good fit of this expansion with our 
numerical results to be presented below in Fig. 6, and 
in view of the interesting bending structure of the nu- 
merically derived current-voltage characteristics in Fig. 
9 below, it seemed promising to calculate the next term 
of the expansion of order 



2 a' {a£Tf 



(53) 



The function /s can also be calculated fully analytically 
and along the same lines: first L2 (x) is derived from (38) 
and (42) and inserted into the o.d.e. (41) for £"2(2;). The 
equation is solved, and the constant of integration is de- 
termined by (45). Then J^^ ^3(0;) dx is derived from (46) 
and the explicit expressions for £i{x) and £2{x) and in- 
serted into (47). However, the result of this calculation is 
still considerably longer than (50) and does not show any 
simple structure. Wc therefore will not give the explicit 
form here. Rather, we summarize essential steps of the 
calculation in more detail in Appendix B. The final steps 
are done by computer algebra (mathematica). Some of 
the results for are shown in Fig. 6 and compared with 
numerical solutions of the full problem (18) - (22). 
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D. Discussion of the result 

Translating back from dimensionless to physical units, 
the result reads 



U{J) = Ut 



+ 



with the dimensionless coefficient 





y Ed%a 








J 2dEa 




{a Erf 


' ) 


' h 


+ . 






{ElApf 





Edla 



2 dEd 



£t a" _ sE^ - (s + 1)E^ 



2 a' 



2E?r 



(54) 



The coefficient of changes sign for £t = Et/Eq = 
[s/(s + With the help of (28), this transition at 



a 



can be located on the Paschen curve; it occurs at 

. l+l/s 



''crit 



(55) 



This is always on the right branch of the Paschen curve, 
since the minimum is at L = e^l^ . 

The result agrees qualitatively with the one given by 
Raizer [21] and Engel and Steenbcck [13]. In particular, 
the leading order correction is also of order a"{j/ jif'. 
However, the explicit coefficient of p differs: while the 
coefficient in [13,21] does not depend on 7 at all, we find 
that the dependence on 7 is essential, as the plot of F in 
Fig. 1 clearly indicates. In fact, within the relevant range 
of 10"^ < 7 < 10°, this coefficient varies by almost four 
orders of magnitude. We remark that it indeed would be 
quite a surprising mathematical result if the Townsend 
limit itself would depend on 7 as in (28), but the small 
current expansion about it would not. Our 7-dependent 
analytical result also excellently fits our numerical solu- 
tions, as we will show in the next section. 



IV. NUMERICAL SOLUTIONS 

We now discuss our numerical results for the voltage u 
as a function of total current j, secondary emission coef- 
ficient 7, mobility ratio n and system size L, as resulting 
from Eqs. (18)-(22). We will work with the Townsend 
approximation a{£) = e"^/'^' (6) with s = 1 as the stan- 
dard case [21]. 



A. The numerical method 

In Appendix A, we showed that the solution u = u{j) 
is unique for fixed 7, /i and L, and we proved the use- 
ful property that the system size L is a monotonically 
decreasing function of the electric field £{0) at the an- 
ode, dL/dSiQ) < (A4), for fixed 7, fi and j. The sec- 
ond observation lays the basis for our numerical iteration 
procedure: 



First the two o.d.e.'s (18), (19) are integrated from 
X = with the known initial value je(0) = j ^-iid some 
guessed initial value £{0) towards larger x. The equa- 
tions are integrated until for some x = x, we find the 
value je{x) = j e~^^ that should be assumed at the fixed 
system size x = L. li x > L, a larger value of £{0) is cho- 
sen for the next iteration step, and if a; < L, a smaller 
£{0) where a linear interpolation of dx/d£{0) is used. 
This iteration loop is continued until the boundary con- 
dition (20) at L is obeyed with sufficient accuracy. The 
potential <p{x) is integrated together with je{x) and £{x) 
by adding the third o.d.e. dx(t> = —£ [34]. The voltage u 
over the system is u = (/>(0) — (^{L). 

For the numerical integration of the o.d.e.'s, we used 
the Isodar.f routine of the ODEPACK package from the 
free-ware site netlib.org. It integrates initial value prob- 
lems for sets of first order o.d.e.'s and chooses automati- 
cally the appropriate numerical method for stiff or non- 
stiff systems. At the same time, it locates the roots 
of any specified function. We defined this function as 
je{x) — j e~^T which returns the value x for the next 
iteration loop with high precision. 



B. Parameters L, 7 and fj,, and j//u-scaling 

The problem depends on the following parameters: the 
first one is the system size L which is proportional to pd 
in physical units. It can take arbitrary values; we explore 
a continuous range of L on both the left and the right 
branch of the Paschen curve. 

The second parameter is the secondary emission coef- 
ficient 7 which is determined by both the gas and the 
cathode surface. Increasing 7 decreases the minimum 
breakdown voltage which is eLj as discussed after Eq. 
(55). This mechanism can be used for improving perfor- 
mance in technical applications like plasma display pan- 
els [33]. According to [21], 7 can take values between 
10~^ and 10~^, in extreme cases even larger. We show 
results either for the two extreme cases 10~^ and 10~^, 
or we show one representative result for 7 = 10~^. 

The third parameter is the mobility ratio /x = /i+//Xe 
of the charged species. Since ions are much heavier 
than electrons, fi is always much smaller than 1. The 
largest value oi fi — 0.0095 [21] is reached for the lightest 
molecules, namely hydrogen. As a standard, we use the 
value /Lt = 0.0035 for nitrogen. 

The functional form of the defining equations (18)-(21) 
and of the small current expansion (53) suggest that u 
in leading order does not depend on j and fj, separately, 
but only on the scaling variable and on the factor 
(1 + /i) w 1. This observation motivates our choice of 
the variable j/ fi in the following figures. However, Fig. 7 
will show that for large j/n'm the abnormal glow regime 
and for large systems L, there is some small /i-dcpendent 
correction to this scaling behavior. Reconsidering (18)- 
(21), this means that the factor (1 -|- fj.) can not simply 
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be equated with 1 even for fi < 10~^, but yields some 
correction. In physical terms, the substitution of (1 + /i) 
by 1 means the elimination of the field increase in the 
anode fall region, and we conclude that in large systems 
in the glow regime, the anode fall yields some small con- 
tribution to the current-voltage-characteristics. 



C. General features of the 
current-voltage-characteristics 



We now give an overview over our numerical results 
in the full parameter regime of the current-voltage- 
characteristics u{j) from Townsend up to abnormal glow 
discharge as a function of rescaled current j/n, system 
size L and secondary emission coefficient 7. In Figs. 2 
and 3, we plot u as a function of j/ fx and L for 7 = 10~^ 
and 7 = 10~^, respectively. The plots follow the style of 
an experimental plot in [25], which is reproduced as Fig. 
1 in [29]. To the best of our knowledge, our Figs. 2 and 
3 for the first time present numerical results in the same 
style. 

Comparing the two figures for different 7, it can be 
noted that on the one hand, the shapes look qualitatively 
similar, while on the other hand, the actual parameter 
regimes of potentials, currents and system sizes vary by 
an order of magnitude or more. Let us now consider the 
common features. 

In the limit of small current j (i.e., in the foreground of 
the figures), the curves saturate to a plateau value which 
actually reproduces the Paschen curve u = ut{L, 7) from 
Eq. (27). 




FIG. 2. M as a function of j/fi and L for the small sec- 
ondary emission coefficient 7 — 10^^. The parameter 
range is 3 • 10"^//i < j/fj. < 5 • 10"^//i for fj, = 0.0035 
and 17.3 <L< 160. 




FIG. 3. Plot as in Fig. 2, but now for 7 = 10"^ 
The parameter range is 10~^ /n < j/n < 7 • 10~^//Lt for 
M = 0.0035 and 3 < L < 28. 

Following u = u{j) along a line of fixed system size L, 
we get the current-voltage-cliaracteristics for this partic- 
ular system characterized by the two parameters L and 
7. For these curves u = u{j), two features can be noted. 

First, for larger L, the voltage u first decreases for 
increasing current j. This is the familiar Townsend- 
to-glow-transition with negative differential conductivity. 
For larger j, the voltage u increases again towards the 
regime of abnormal glow. However, in the minimum of 
the potential, there is no plateau in contrast to experi- 
mental plots. This is because we solve the purely one- 
dimensional system without the possibility of a lateral 
growth of the glow discharge column. 

Second, for smaller values of L, in particular, when 
starting from the left branch of the Paschen curve, the 
voltage does not decrease for increasing current, but it 
increases immediately. We will discuss this different bi- 
furcation structure in more detail at the end of this sec- 
tion. 



D. Spatial profiles 

It is instructive to study the spatial profiles of electron 
current je{x) and field £{x) for different system sizes. 
In Figs. 4 and 5, we plot such profiles for L — eL^ 
and for L = e^L^ = eLcrit- The smaller system size 
eL^^ coincides with the minimum of the Paschen curve, 
while Lcrit = e^L-y (55) is the system size where a" 
in (53) changes sign. So for L < Lcrit the voltage in- 
creases initially in the small current expansion around 
the Townsend limit, while for L > Lcrit it decreases. 

Note that in contrast to previous plots, e.g., in [21], our 
cathode is on the right hand site at x — L, because we 
found it more convenient to work with a positive field £. 
The electron current is normalized by the total current. 

In each plot, the profiles for the two smallest cur- 
rent values are well described by the small current ex- 
pansion from Section III. This is in agreement with the 
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range of validity of these expansions of < 0.08 for 
L = eL^ = Lcrit/e or of j/jj, < 1.3 • 10"^ for L = eLcriu 
resp., estimated according to (52). 




Q r^^^^ ■ ^— ■ ' ' ' >- 

2 4 6 8 10 „ 12 



FIG. 4. Spatial profiles je{x)/j and £{x) for system size 
L = eL~f at the minimum of the Paschen curve. Plotted 
are curves for j/ n = 0.01, 0.1, 0.3, 1, 3. Other parame- 
ters: 7 = 0.01, n = 0.0035. 
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FIG. 5. The same as in the previous figure, but now for 
the larger system size L = e^L^. The current now ex- 
plores the smaller values j/fi = 10"", 10"^, 3-10"^, 0.01, 
0.03, 0.1, 0.3. 

For larger currents, a separation into resistive column 
on the left and cathode fall on the right becomes pro- 
nounced. While in the smaller system, both regions take 
about equal parts, in the larger system, the cathode fall 
takes only a small part of the volume on the right hand 
side. 



E. Comparison of numerical and analytical results 

Let us now compare the current- voltage-characteristics 
corresponding with these profiles with our analytical re- 
sults from Section III. In Fig. 6, the numerical results for 
u{j) are plotted as a thick solid line, and the analytical 
expansions (51) and (53) up to second or third order in 
j/fj, as thin solid and dashed lines, respectively. For the 
calculation of the third order, the procedure described 
in Appendix A has been followed. Fig. 6 shows that in 
particular the expansion up to order (j//Lt)^ gives a very 
good agreement at least within the range of validity of 
j/fj, < 0.08 or 1.3 • 10~^, resp., according to (52). 




FIG. 6. u{j) for the systems from Figs. 4 and 5 in the 
small current limit: upper plot L = eL-y, lower plot 
L = e^L-y, both with 7 = 0.01. Numerical result (thick 
solid lino), analytical result (51) up to second order in 
j/n (thin solid line) and analytical result (53) up to third 
order in j/fj, (dashed line). 



F. Corrections to j//U-scaling 

Fig. 7 shows the currcnt-voltagc-characteristics for the 
same two systems, but now up to larger values of the 
current than in Fig. 6. Actually, the same current range 
is explored in each system as in the corresponding Figs. 
4 and 5. 

In addition, in Fig. 7 we test the j/ /^-scaling by plot- 
ting w as a function of j7/x within the physical range of 
/i-vahics, including the limit of /i = 0. It can be noted 
that for short systems or small currents, the /x-correction 
is negligible; this means that (1 -|- /x) can be replaced by 
1 in Eq. (19) without visible consequences. In contrast, 
for large systems and large currents, there is a small, but 
visible /U-correction to the dominant j/ /x-scaling. 
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25 




FIG. 7. Current-voltage characteristics for the same two 

systems as in Figs. 4, 5, and 6 (upper plot L = eL~f, 
lower plot L — e^L^), but now for a larger current range 
than in Fig. 6. Furthermore, besides the curves for the 
nitrogen value for the mobility ratio /i = 0.0035, also the 
curves for fj, = 0.0095 (hydrogen), n = 0.001 and the 
limiting value /^ = are shown. 



G. Discussion of bifurcation structures 

Wc now set the final step in tlie quantitative un- 
derstanding of the current-voltage-characteristics at the 
transition from Townsend to glow discharge. We charac- 
terize the transitions as subcritical, mixed or supercritical 
and locate them in parameter space. 

Fig. 8 gives an overview over the different behaviors for 
7 = 0.01. It corresponds to different iy-sections of plots 
as in Figs. 2 and 3, but now with plotted on a linear 
rather than a logarithmic scale. For the terminology of 
sub- or supercritical bifurcations, it should be noted that 
= is a solution for arbitrary u. So the complete 
u-axis is a solution, too. 

In the case of L = 0.85Lcrit) there is a pure forward 
or supercritical bifurcation: u increases monotonically as 
increases. In contrast, for L = l.OBLcrit, the bifurca- 
tion is purely subcritical: with increasing the voltage 
first decreases, and eventually it increases again. This 
subcritical behavior continues down to L = Lc-iL = L^e^ 
where a" changes sign. So indeed, the sign change of a" 
in the small current expansion determines the transition 
from subcritical to some other behavior. However, for 
L < Lcrit, the system not immediately enters the super- 
critical regime, but some imexpcctcd "mixed" behavior 
can be seen: for increasing j//x, the voltage u first in- 
creases, then it decreases and then it increases again. We 
distinguish "mixi" where the voltage minimum at finite 
j//Lt is smaller than the Townsend voltage, and "mixn" 
where it is larger. 




0.05 0.1 0.15 0.2 . , 0.25 

FIG. 8. The current-voltage-characteristics for fixed pa- 
rameters 7 = 0.01 and n = 0.0035 and different system 
sizes, measured in multiples of Lcrit = L^e^. Shown are 
all possible bifurcation structures from supercritical up 
to the familiar subcritical case for various values of L. 



Fig. 9 shows a zoom into Fig. 8: a smaller range of 
current and of system sizes. The form of an upwards 
parabola of the order {j/jJ,}^ next to the w-axis is well 
described by the analytical small current expansion of 
Section III. However, in contrast to initial hopes, the 
turn-over of the curves at j/jj, 0.02 is not covered by 
the small current expansion up to order {j / ^if whose 
coefficient even changes sign within the parameter range 
of Fig. 9. In fact, the range of validity of the expansion 
breaks down at < 0.017, just briefly before the first 
interesting bending structure in the characteristics. 




0.02 0.04 0.06 0.08 i / n 0.1 



FIG. 9. Zoom into Fig. 8 with smaller values of j / ji and 
system sizes L. The range of validity of the analytical 
small current expansion (53) is j/ n < 0.017. 



10 



Finally, in Fig. 10, the bifurcation behavior in the full 
parameter range of 7 is explored. The transition from 
subcritical to mixj always takes place when a" changes 



sign, i.e., at system size 



L^e . The transitions 



from mixi to mixn and then further to supercritical oc- 
cur at smaller relative system sizes L/Lcrit when the sec- 
ondary emission coefficient 7 is smaller. All transitions 
occur on the right branch of the Paschen curve, since its 
minimum is at L/Lcrit = = 0.368. 
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FIG. 10. Complete overview over the bifurcation behav- 
ior of the current-voltage characteristics £is a function of 
secondary emission coefficient 7 and system size L. 



V. SUMMARY AND OUTLOOK 



[13], but with a different, strongly 7-dependent propor- 
tionality constant. 'We also have calculated the term of 
order These analytical expansions are in very 

good agreement with our numerical results within their 
predicted range of validity. 

Of course, the study of this minimal model can only be 
a first step, and it has been suggested to include a num- 
ber of additional features. First, 7 might not be constant; 
while the dependence 7 = 7(1, V) [18] on global parame- 
ters seems unphysical, the local dependence 7 = ^{E/p) 
has been suggested [35] and experimentally tested [36]. 
Second, the particle mobilities might be field-dependent 
IJL± = n±{E). Third, diffusion was neglected, the ap- 
proach was fully local, and only one ion type was consid- 
ered. However, our aim was first to settle the predictions 
of the classical model in full parameter space as a corner 
stone and starting point for any future extension that 
simultaneously also will increase the number of parame- 
ters. 

The motivation for this work is the impressive vari- 
ety of spatio-temporal patterns formed in short barrier 
discharges [4 12]. The nonlinear element responsible for 
the spontaneous pattern formation is believed to be a 
gas discharge in the parameter range of the present work. 
Parameter regions with negative differential conductivity 
(NDC) are generally believed to play a decisive role in the 
formation of the instabilities. Knowledge about NDC re- 
gions and the bifurcation structure in the range of these 
experiments therefore are a condition for their future 
investigation, and conversely, properties of the current- 
voltage characteristics might be deduced from temporal 
oscillations or current constrictions. These pattern for- 
mation processes will be subject of our future studies. 



We have studied the classical minimal model that cre- 
ates a Townsend or glow discharge, in one-dimensional 
approximation. The dimensionless current-voltage char- 
acteristics u = u{j/fi) depends on essentially only two 
parameters, the secondary emission coefficient 7 and the 
dimensionless system size L oc pd. ('V\^itli j/^-scaling, 
the further dependence on the small mobility ratio fi = 
^i+IHe is very weak and becomes visible only for long 
systems in the normal and abnormal glow regime.) Nu- 
merically, we have fully explored the bifurcation structure 
as a function of 7 and L. Besides the familiar subcriti- 
cal bifurcation structure of long systems, for decreasing 
system size L, there is a sequence of current- voltage char- 
acteristics, that we have called mixj and mixn, before the 
supercritical transition is reached. This general sequence 
is the same for all relevant values of 7 while the precise 
lengths where the transitions occur, depend on 7, cf. Fig. 
10. Analytically, we have calculated the small current 
expansion about the Townsend limit in a systematic ex- 
pansion. 'We found that the term of order j/ jj, is missing, 
and that the term of order (j /fJ.)'^ indeed is proportional 
to the second derivative of the Townsend coefficient a" 
according to the old argument from Engel and Steenbeck 
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APPENDIX A: UNIQUENESS OF THE 
SOLUTION OF THE BOUNDARY VALUE 
PROBLEM 

Here we prove that the boundary value problem de- 
fined by (18) - (21) for fixed j, /z, 7 and L defines a 
unique solution (7e(.x), £(.t)) and hence a unique poten- 
tial u = Jq £{x)dx. This lays the ground for our analyt- 
ical as well as for our numerical procedure. We will keep 
J, /i and 7 fixed within this appendix, and will discuss 
how £{x) is determined by L and vice versa. 

First Eq. (18) for je is integrated with initial condition 
jg(0) =j (20) and inserted into (19): 

^£d^£ = j^-{l + ^) r «(^(-)) (Al) 
The boundary condition (20) at L amounts to 
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a{£{x)) dx = , 



(A2) £2ix) 



where we assume that 

a(f ) > and da/d£ > for all f > . (A3) 

An initial condition £{0) defines a unique solution £{x) 
of (Al) and hence a unique system size L through (A2). 
We will show below that i is a monotonically decreasing 
function of £{0) 



dL/d£{0) < for fixed j, ji, 7 . 



(A4) 



This statement has two immediate consequences: (i) 
it shows that £{Q) and therefore also £{x) and u are 
uniquely determined by L; and {ii) it lays the ground for 
our numerical iteration procedure where £"(0) is fixed, 
and the resulting L is calculated and compared to the 
true L. 

Why is statement (A4) true? Compare two solutions 
£i,2{x) and suppose that £i(x') > £2(2;') on some inter- 
val < x' < X. Then for the difl[erence, we get from (Al) 
that 



{d,£^ - d,£l) 



2j (1 + ^) 

_ g- "(£2(3;)) dx _ ^- a(fi(x)) da; ^ Q 



(A5) 



where the bound . . . > is a direct consequence of (A3). 
So when £1 is above £2 on some interval < x' < x, then 
at the end of the interval, we have dx£i > 9x^1 > ^^id £1 
stays above £2- As a consequence 

£i{x) > £2{x) for all a; > , if ^:i(0) > £2(0) . (A6) 



2a2^2£:3 



(1 + m) 



o! £t 



- {axye-""' + (1 + At)(e-^""' - 26""^) 



2(':^-l + (l + M)^ 



(aa; + l)e-"^| 



—ax [ Lj — ax + 2(1 + /x) — 

+2(1 + M) (^^ + (1 + ^)1^ 



ax e 



-(1 + M)'e 



-2ax 



+ c 



(B2) 



The constant of integration C is determined through (45) 



/' 

Jo 



£2 dx 



2 a' Q3^i2£2 



(B3) 



Finally, /q £3 dx is determined by £1 and £2 through 

(46) as 

/ £3dx = -— £i£2dx--— £i dx. (B4) 
Jo a Jo 3! a' Jo 

Insertion of the result in Eq. (47) yields the third order 
expansion (53) with an explicit expression for the func- 
tion /s. 

All integrals can be performed analytically. However, 
the results are lengthy and exhibit no simplifying struc- 
ture. Therefore, we rather have performed the remain- 
ing integrals by computer algebra. Results are shown in 
Fig. 6. 



Inserting this into (A2) and using (A3), statement (A4) 
results. 



APPENDIX B: THE CORRECTION OF 0{J^) 
ABOUT THE TOWNSEND LIMIT 



We here sketch the essential elements of the calculation 
of the third term of the expansion about the Townsend 
limit: first Eq. (38) is integrated. With (48) for li, (49) 
for £1, and the boundary condition (42), we get 



l'2{x) 



a'^^£j 



{axY 



— ax 



(l + M)(l-e-«-) 



(Bl) 



This result allows us now to integrate £2(2;) in (41) with 
the rather lengthy result 
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